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ABSTRACT 

We study the limits of accuracy for weak lensing maps of dark matter using dif- 
fuse 21-cm radiation from the pre-reionization epoch using simulations. We improve 
on previous "optimal" quadratic lensing estimators by using shear and convergence 
instead of deflection angles. We find that non-Gaussianity provides a limit to the accu- 
racy of weak lensing reconstruction, even if instrumental noise is reduced to zero. The 
best reconstruction result is equivalent to Gaussian sources with effective independent 
cell of side length 2.0/i^^ Mpc. Using a source full map from z=10-20, this limiting 
sensitivity allows mapping of dark matter at a Signal-to-Noise ratio (S/N) greater 
than 1 out to I < 6000, which is better than any other proposed technique for large 
area weak lensing mapping. 

Key words: Cosmology-theory-simulation-observation: gravitational lensing, dark 
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1 INTRODUCTION 

The lens mapping of dark matter is an essential cornerstone of modern precision cosmology. Weak gravitational lensing has 
developed rapidly over the past years, which allows the measurement of the projected dark matter density along arbitrary 
lines-of-sight using galaxies as sources. Recently, Smith et al. (2007) have demonstrated the first CMB lensing detection. The 
goal is now to achieve high precision cosmological measurements through lensing, at better than 1% accuracy. 

Galaxies are plentiful on the sky, but their intrinsic properties are not understood from first principles, and must be 
measured from the data. Future surveys may map as many as 10^" source objects. Using galaxies as lensing sources has several 
potential limits (Hirata & Seljak 2004), including the need to calibrate redshift space distributions and PSF corrections, to 
be better than the desired accuracy, say 1%. This will be challenging for the next generation of experiments. 

Some sources, such as the CMB, are in principle very clean, since its redshift and statistical properties are well understood. 
Unfortunately, there is only one 2-D CMB sky with an exponential damping at I ^ 1000, which limits the number of source 
modes to ~ 10^. 

The potential of detecting the 21-cm background from the dark ages will open a new window for cosmological detections. 
Studying the 21-cm background as high redshifts lensing source, as well as the physics of the 21-cm background itself, provide 
rich and valuable information to the evolution of universe. The number of modes on the sky is potentially very large, with 
numbers of 10^^ or more. For this reason, 21-cm lensing has recently attracted attention. However, most of the reconstruction 
methods are based on a Gaussian assumption (Pen 2004; Cooray 2004; Zahn & Zaldarriaga 2006; Benton Metcalf & White 
2006; Hilbert et al. 2007). In contrast to CMB lensing, where the Gaussian assumption works well, non-Gaussianity in 21-cm 
lensing may affect the results. Non-linear gravitational clustering leads to non-Gaussianity, and ultimately to reionization. In 
this paper, we will address the problem of the lensing of pre-reionization gas. 

21-cm emission is similar to CMB: both are diffuse backgrounds. It is natural to apply the techniques used in CMB 
lensing. Hu & Okamoto (2002) expand the CMB lensing field in terms of the gravitational potential (or deflection angles), 
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and construct a trispectrum based quadratic estimator of potential with maximum S/N. However, unlike CMB, the 21-cm 
background has a 3-D distribution and is intrinsically non-Gaussian. A fully 3-D analysis is explored in Zahn & Zaldarriaga 
(2006), where they generalize the 2-D quadratic estimator of CMB lensing (Hu & Okamoto 2002) to the 3-D Optimal Quadratic 
Deflection Estimator (OQDE). 

A local estimator was proposed in Pen (2004), which assumed a power law density power spectrum. In this paper, we will 
design localized estimators for the lensing fields under the Gaussian assumption, and apply the derived reconstruction technique 
to Gaussian and non-Gaussian sources. The influence of non-Gaussianity can be measured by comparing the numerical results 
between the Gaussian sources and non-Gaussian sources. 

Quadratic lensing reconstruction is a two point function of the Icnsed brightness temperature field of the 21-cm emission. 
In the paper, 3-D quadratic estimators are constructed for the convergence (k), as well as the shear (7). Our method recovers 
the K and 7 directly instead of gravitational potential or deflection angles. Our estimators have in principle the same form as 
the OQDE, consisting of the covariance of two filtered temperature maps. The OQDE reconstructs the deflection angle, while 
our estimators reconstruct the kappa and shear fields. Our filtering process can be written as a convolution of the observed 
fields. As presented in Appendix and section 4, our combined estimator is unbiased, and equally optimal as the OQDE for 
Gaussian sources, and has better performance for non-Gaussian sources, and recovers three extra (constant) modes. 

Other authors also developed reconstruction methods from alternative approaches. Benton Metcalf & White (2006) give 
a estimator for shear. They choose the separate 2-D slices at certain redshift intervals, and then these slices can be treated as 
independent samples for the same lensing structure. As a result, the information between these slices are lost. Cooray (2004) 
expands the lensed field to a higher order of the gravitational potential, and investigates the higher order correction to the 
lensed power spectrum. 

The paper is organized as follows: The basic framework of lensing and the reconstruction method is introduced in §2. 
The numerical methods are presented in §3. The results are discussed in §4. We conclude in §5. 



2 LENSING AND RECONSTRUCTION 

Photons are deflected by clumpy matter when they propagate from the source to the observer. This effect can be used to map 
the mass distribution if we can measure the distortion of an image. In this section, we will first review the lensing theory, 
which serves to define our notation. We then develop an optimal quadratic estimator using a maximum likelihood method. 
The reconstruction depends on the power spectrum of the source. The noise and normalization of the reconstruction are 
calculated in the appendix. 



2.1 Lensing 

The Jacobian matrix describing the mapping between the source and image planes is defined as 

Here x is the radial coordinate, and /k (x) is the comoving angular diameter distance. We consider a ray bundle intersecting 

at the observer and denote x{G,x) the comoving transverse coordinate of a ray. 

In the lensing literature, the physical quantities frequently used to describe a lensing field are convergence k and shear 
7, which are given by 



J(0,x) 



1 — K — 71 —72 

-72 1 - /t -h 71 



Equivalently, the convergence and shear can also be written as k = ($,11 + $,22)/2 ; 71 = ($,11 — $,22)/2 ; 72 = $,12- ^ is 
the projected 2-D potential: 



* = - 



Jo JK(X) 

Here subscripts '1' and '2' refer to the derivative to the two perpendicular transverse coordinates, and is the 3-D Newtonian 
gravitational potential. Note that the integral is along the actual perturbed path of each photon. In the Born approximation, 
the defiection is approximated by an integral along the unperturbed path. 

In the small angle approximation (Limber 1954), Vj_ can be replaced by in the integral. We get the Limber equation 



f dxW,x)^, (3) 
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with g{x' , x) = fK{x')fK{x ~ x')//^(x)- is the mass density parameter, Ho is the current Hubble constant, a is the scale 
factor, and S is the over-density. 

Kaiser (1992) derived the Fourier-space version of the Limber equation 

Here Pk{1) is the 2-D power spectrum of the k field, P{1/ fK{x)^x) is the 3-D power spectrum of matter, and xh is the comoving 
distance to the Hubble horizon. The equation is valid when the power spectrum Pk evolves slowly over time corresponding to 
the scales of fluctuation of interest, and these fluctuation scales are smaller than the horizon scale. 



2.2 Reconstruction of large-scale structure 

We first heuristically review the quadratic lensing estimation in two dimensions. Then we will proceed with a generalization 
to 3-D with a quantitative derivation. 

Lensing will change the distribution of a temperature field by changing length scales. Lensing estimation relics on statistical 
changes to quadratic quantities in the source plane temperature field. We use a tilde to denote a lensed quantity. All estimators 
work by convolving the temperature field with a window, 

fi(x) = I d^xf{x')Wi{x-x') , (5) 

and a second window 

f2(x) = J d^x'f(x')W2{x - x') . (6) 
The quadratic estimator is simply the product of the two convolved temperature fields, 

E{x) = fi{x)f2ix). (7) 

In the weak lensing case, the estimator is a linear function of the weak lensing parameters (k, 7). The simplest case is two 
equal, azimuthally symmetric window functions Wi = W2 = fir). Considering the limit that k is a constant value, the 
estimator is linearly proportionate to k: 

{E}cxK + V , (8) 

and V is the mean covariance. Here (...) means ensemble average. For a stochastic random field, the ensemble average can be 
calculated by the volume average if the volume is big enough. We can absorb V as well as the normalization coefficient into 
E for convenience, i.e., E{x) = Ti{x)T2{x) — V . When k is spatially variable, E needs to be normalized by a scale dependent 
factor h{k). This corresponds to a convolution of n with a kernel: 

{E{x)) = J d''x'K.{x')b{x - x') , (9) 

where kernel b is the Fourier transform of the normalization factor. 

One can optimize the functions to minimize the error ou the Icnsiug variables. In this paper we will compare various 
forms of the smoothing windows, which include as special case the traditional Optimal Quadratic Deflection Estimator. The 
simplest case is a constant value of k, for which one can compute its variance 

(k^) = mixff). (10) 

Lensing is a small perturbation of the variance, therefore we can calculate the variance from the unlensed source field, i.e., 
{{Ti(xff) ^ {{fi{xff). Performing a variation to minimize the variance, one can find the optimal window function. It 
turns out that the window functions do not depend on the spatial structure of the lensing field. Only the normalization factor 
b in Eq. (9) is scale dependent. We solve the optimal window function at scales where the constant k approximation works 
well, and the solution should also be optimal for other scales. 

Shear and deflection angles are tensorial and vectorial quantities and require anisotropic or vectorial choices of the window 
function 

= T1T2, Ti = J d'^9'f{e')Wi{0 - 6'), T2 = J d'9'f{e')W2i0 - 6') , (11) 

Ed = Tifa, Ti = j d^e'f{e')w-i.{0 - e'), fa = j A'e'f{e')W2{e - e') . (12) 

This will be explained in detail in sections 2.2.2 and 2.3. 

The source is usually treated as a Gaussian stochastic field in the literature on reconstruction methods. While this is valid 
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for CMB on large angular scales, 21-cm background sources are not always Gaussian. In this paper we attempt to understand 
the influence of this non-Gaussianity. Optimal estimators for Gaussian sources are not necessarily optimal for non-Gaussian 
sources. Here, wo will construct the convergence and shear field directly, instead of following the deflection angles or potential 
fleld reconstruction in CMB lensing. There are three reasons to do this: Firstly, the strength of lensing is evident through the 
magnitude of k or 7 since they are dimensionless quantities. The rms deflection angle of photons from 21-cm emission is at the 
magnitude of a few arcmin, which is comparable to the lensing scales we arc resolving. Some authors argued that perturbation 
theory on the deflection angle will break down at these scales (Cooray 2004; Mandel & Zaldarriaga 2006). However, k and 7 
are still small and can still work with perturbation calculations without ambiguity. Secondly, k and 7 have well deflned limits 
as they approach a constant, while only spatially variable deflection angles or potentials can be measured. This signiflcantly 
simplifles the derivations. Finally, k and 7 are standard variables to use in broader lensing studies, such as strong lensing and 
cosmic shear. Using the same convention in different subfields will help to generalize the underlying physics of lensing. 

The estimators are unbiased, as shown in the appendix. Furthermore, we conflrm that our combined estimators from k 
and 7 have the same optimality as the OQDE for Gaussian sources. When the sources are non-Gaussian, our estimators have 
better S/N. 



2.2.1 Maximum likelihood estimator of k 

We now derive the quantitative window functions for 21-cm lensing reconstruction. Due to their similarity, it is helpful 
to quickly review the reconstruction in CMB lensing: The early work by Zaldarriaga & Seljak (1999) used the quadratic 
combination of the derivatives of the CMB field to reconstruct the lens distribution. Since the CMB has an intrinsic Gaussian 
distribution, The optimal quadratic estimator (Hu 2001b) can also be applied to lensing reconstruction with CMB polarization 
(Hu & Okamoto 2002). Zahn & Zaldarriaga (2006) generalized the optimal quadratic estimator of CMB lensing to 21-cm 
lensing. 

We will construct estimators for k and 7 with the 21-cm brightness temperature fields, starting from a maximum likelihood 

method which is consistent with the quadratic minimum variance method when the field is Gaussian. We will show that the 
OQDE and our approach are the same if the sources are Gaussian, however the problem is simplified in a intuitive way by 
using the limit that k and 7 vary slowly in small scales. 
The magnification is 

^' = n V2 2 ~ 1 + 2« . (13) 

(1 — k)^ — 7 

The last approximation is valid since both k and 7 are much smaller than 1 in the weak lensing regime. 

We use Bayosian statistics and assume the prior distribution of parameter k to be flat. For a M pixel map on the sky, 
the posterior likelihood function of the source field has a Gaussian distribution, and can be written as 

r{f(k)) = {27v)-^'^''dct{Cff)'h-i'^^'^ff''^ . (14) 

Here T = Tj, + n is the brightness temperature of the difiiusive 21-cm emission lensed by the large-scale structure plus 
measurement noise. To simplify the algebra, we use the negative logarithm jC of the likelihood function in our calculation, 

jC = -\nP = ^f^Cflf + ^\ndetCff. (15) 

Here T is the 3-D discrete Fourier transform of measured temperature. C^.^. = Cs + Cn is the covariance matrix, and the 
signal contribution Cs and noise contribution Cn arc both diagonal in Fourier space and uncorrelated to each other. In the 
continuum limit, the likelihood function can be written as 



d^fcEi^j . (16) 



We use 7^*d' = P-juik) + Pnik) to represent signal plus noise power spectrum in the following text, where Psoik) is 3-D power 
spectrum of the distorted 21-cm field, and P^ik) is the noise power spectrum. The dimensionless power spectrum of the 3-D 
21-cm gas slices can be written as 

^3D(fc) = ^PMk) , (17) 

where k = \k\ since the gas is statistically isotropic. 

The geometry of the 21-cm field will be changed by the lensing: 

Tb(fcx,fc||) = j d\n{x)e-"' '' = Jd\± y"dx||r(,((l-K)a;x,X||)e-'('=^ "^+*il"ll' = ■^^-i^T,((l + «)fcx,fc||) , (18) 
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where '_L' and '||' mean the perpendicular and parallel direction of the line-of-sight respectively. We ignore the contribution 
of shear first. Then the length scale is magnified on the transverse plane by a factor k. Isotropy is broken in 3-D but is still 
conserved on the 2-D cross section. The statistical properties of the 21-cm field will be changed by the lensing, i.e., the power 
spectrum will also change: 

{f;(fcx,fc|t)T'6(fcY,fc;|)) = (2^)=^5"°(fcx-fe'x)(27r)<5°(fc|| -fcf|)P3D(fex,fci|) . (19) 
The delta function has the property: 

+ K)fex - (1 + K)fel) = 7T^^5'°(fex - fel) . (20) 

(1 + K) 

Therefore the relationship between the unlensed and lensed power spectrum is 

P3D(fex,A;||) = (l-|-2K)P3D((l-|-K)fex,A;||) = (1 + 2k)P3d(^(1 + K)^kl + W (1 + 2/t)(P3D(fc) + /tAP3D) , (21) 

where AP3D = Psok^k^/k'^), and P^oik) = dP3D{k)/dk. The second equivalence is due to the statistical isotropy of the 

unlensed power spectrum. 

Differentiation of the lensed power spectrum gives 

= 2P3D + (1 + 2k) AP3D , (22) 
and the maximum likelihood condition requires 
SC 1,3 /• d^fc (Pjg - |f pL-3) ^P3D 

Sk~ J (2^)3 ptct2 Sk ^ 

which has the solution 



/ 



(277)3 



{\ffL-^)J'^{k) - K . (24) 



We have approximated P^$ by P30' in the denominator of Eq. (23) . To simplify the problem, we assume the source is a cube 
with physical length L in each dimension. The offset constant V„ = {a^) = J k / {2n)^ P^^ {k)J^'^ {k) , and the optimal filter 
JF-^ is 

2P3 D(fc)+AP3D (fc) 

^3D* {k)QK 

with Q« = /d»fc/(27r)''(2P3D + AP3D)(fc)J^"(fc). 

From Parseval's theorem, we can rewrite Eq. (24) in the form of a convolution of the density field and a window function 
in real space 



I 



d-^k 
(2^ 



f*{k)f{k)J^''{k) = J d'xf:^{x)f^^{x)=L^ J da:||f,:^(xx,a;|,)r^,(xx,a;||). (26) 



In Eq. (26) the two window functions are the decomposition of the optimal filter Wi{k)W2{k) = J^"'{k). The last '=' in 
Eq. (26) holds when /t is constant. One can choose W^{k) = W5{k) = V^. If < 0, we choose Wf = -W^ = 
The convergence field is equivalent to the covariance of the measured maps with two windows applied. In the slowly spatially 
varying ft limit, all decomposition into two windows are equivalent. As we will show later, the shear construction can also 
be represented in the form of the covariance of two filtered temperature maps. These maps will have symmetric Probability 
Density Function (PDF), which can reduce the non-Gaussianity of the maps so that a better S/N level can be achieved, when 
the shear window functions are chosen properly. The last two steps in Eq. (26) assumes the fluctuation of the convergence 
field is slow compared to the filter. Then we can apply the estimator to each beam in the map: 

£;«(xx) = j dx||f^^(cc)f^Jx) - 14 , (27) 

where T^^ and T^^ are the convolution of T and window function Wi{x) and W2{x) respectively, which are the real space 
version of Wf (fc) and W2{k). The reconstruction of the k is dominated by the gradient of the power spectrum dlnA^/dlnfc, 
which follows the expression of our estimator in Eq. (24). 

We can then generalize the estimator to a spatially varying lensing field. In the appendix we show 

j d'a;'x/t(a;'x)6.(a:x - x'^_) = (E^xi.)) . (28) 

Equivalently, for smaller scales, we will need to normalize the reconstructed lensing field by a scale dependent factor in Fourier 
space, which is calculated in the appendix. 
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k(l) = h-^(l)E4l) = k(0 + n{l) , (29) 

where I = k±x{zs), and Zs is the redshift of the source. Here 6^(0 is the normahzation factor (hmj^o &k(0 = l)i '^(0 is 
the noise, since different Fourier modes are independent. They do not depend on direction because variables related to k are 
isotropic on the transverse plane. In the appendix, wc show that the normalization factor is unity at small I when has the 
form as Q„ = / d^/c/(27r)=^(2P3D + AP3D)(fc).;^"(fc). 



2.2.2 Estimator of shear 

When shear is taken into account, not only the scale but the directions of the coordinates are changed. We will start the 
derivation from the constant shear case. 

T6(fex,/C||) = J d^xft,{x)e-"' '' = y"dVy"da;||r,,(Jxx,X||)e-'('=^-"^+'=ll"ll) 

= \J\-' J d\; J da;||T,(xl,a;||)e-'('=^-^'^+'=ii^i|) = |jr ^T,(J-^fcx, fey) , (30) 

here d"^!'^ — \ J\d'^x±,k'j_ = J~^fcx. Now the symmetry is broken even on the transverse plane due to the anisotropic distortion 
caused by the shear. 

Since <52°(J-ifc) = |J|5^°(fe), Eq. (19) implies 

f'3D(fcx, fc||) = |jrip3D(J"'fcx, fell) w (1 + 2K)[P3D(fe) + AP3o{k){K + 7i cos26>fe^ + 72 sin26ifc^)] , (31) 

where 9kj^ is the angle between fcx and the transverse coordinate. 

Maximum likelihood requires 5C/5^\ = 0, and 5C/5^2 = 0. The maximum likelihood shear estimators can be written as 

a tensor E-> : 



(32) 



where TJ. is convolution of the temperature field with , and {k) = {2APzY>/P^Y)Q-y)^^^ki, fc (i,j = 1,2) is one 
of the two unit vectors on the transverse plane. When AP < 0, we can choose = \2/\Psii/PhQ^,\'-'^ki,W:^ = 

-|2AP3D/P3DQ7r^^fe2. The normalization factor Q-y = / d^fc/(27r)^AP3D(fc)feifc2W^7(fe)lV2'(fe)- The two components of shear 
are now: 

A - P - P A - -^711 ~ -^722 /OQ\ 
7l — tj~)\2 — 1^-)21, 72 — ^ . \o6) 

Note that there is a difference between the reconstruction for convergence and shear. Shear reconstruction depends on 
the gradient of P{k), while convergence reconstruction depends on the gradient of A^(fe) in a 2-D analogue. To tost our 
method, we can generated a Gaussian source field with power law power spectrum P(fe) = fe'^. In the 2-D analogue case, 
the convergence field can not be measured if /3 = —2, because the variance is conserved. However in 3-D, when (3 = —3, the 
convergence field can still be measured, which is due to the more complicated shape of the window function in 3-D. When 
/3 = 0, the shear can not be measured in either 2-D or 3-D. 

In analogy to k reconstruction, we can calculate the normalization factors bj^ and 6^,3 • The calculations for the normal- 
ization factors and noise are presented in the appendix. 



2.3 The combined estimator and the OQDE 

The combined estimator of k can be written as 

/tcombitj- iv4z)-i + iv^jz)-i ' ^•^^^ 

where 7e is the convergence constructed from shear field, 

7e(Z) =7i(Ocos26»i-h72(Z)sin20, , (35) 
and 01 is the angle of I. 

The 2-D OQDE in CMB lensing can be written as product of two filtered temperature field (Hu 2001b; Lewis & Challi- 
nor 2006). Furthermore, the 3-D OQDE can be written in the same form of Eq. (12), though it is not explicit (private 
communication with Oliver Zahn). 



E. 



i{0) = L-^ j dx||Ti(0,X||)r2(0,x,|) , (36) 
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and 

J d^e'd{0')ba{e - e') = {Ea{0)) . (37) 

6d is a normalization factor, Ti = / d^9'T{0')Wi{e - 0') and T2 = / A^e'T{e')W2{e - 6') are convolved temperature fields, 
where the window functions are Fourier transforms of: 
-ilPzT,{l,k\\) 



Wi(f,fe||) 



W2{l,k») = i . (38) 

We note that the OQDE and our estimators have the same form. The contribution from lensing in Eq. (10) is secondary, and 
the noise of reconstruction is mainly determined by the unlensed terms. Therefore we can measure the numerical reconstruction 
noise without lensing the sources. 



3 NUMERICAL METHODS 
3.1 Simulation 

The fluctuation in the 21-cm brightness temperature may depend on many factors, such as the gas density, temperature, 
neutral fraction, radial velocity gradient and Lya flux (Barkana & Loeb 2005). In our work, we do not consider the redshift 
space distortion effect caused by the non-zero radial peculiar velocity gradient, and simply assume the brightness temperature 
is proportional to the density of the neutral gas. 

T,«(27mK)(i±£)'''^i^^(l + &i) , (39) 

where is the brightness temperature increment respective to CMB, is the spin temperature, which is about to be much 

bigger than Tcmb , and fei is the over-density of the neutral hydrogen. 

Our work mainly focuses on the non-Gaussian aspect and 3-D properties of the reconstruction, and these effects also 
exist in a pure dark matter distribution. The neutral gas will trace the total mass distribution, which is dominated by the 
dark matter haloes. A simplification is to use the dark matter as the source directly. Even though this will bring some bias at 
small scales, the approximation is valid at large scales. The dark matter distributions are generated using the PMFAST code 
(Merz et al. 2005). 

The high resolution PMFAST simulation was performed on a 1456^ fine mesh with 3.9 x 10*^ particles. The production 
platform was the IA-64 'lobster' cluster at CITA, which consists of 8 nodes. One of them was upgraded, so we have used the 
remaining 7 nodes. Each node contains four 733 MHz Itanium-1 processors and 64 GB RAM. The simulation started at an 
initial redshift Zi = 100 and ran for 63 steps with comoving box-size L — 50ft~^ Mpc. The initial condition was generated using 
the Zeldovich approximation, and the matter transfer function was calculated using CMBFAST (Seljak & Zaldarriaga 1996). 
The cosmological parameters were chosen in accordance with the WMAP result (Spergel 2003): f2m = 0.27, = 0.73, fib = 
0.044, n = 1.0, (T8 = 0.84, and /lo — 0.71. 20 independent boxes were generated. We had 3-D data at = 7 at hand, and used 
them in our numerical tests for convenience. 



3.2 Convergence and shear map construction 

The dimensionless power spectrum, which is the contribution to the variance of over-density per logarithmic interval in spatial 

wave number, can be measured from the source data in the periodic simulation box. 

To reduce the computation time, our numeric results on the reconstruction use a re-sampled distribution. We will generate 
20 independent sources each on 512^ grids, to investigate the statistics. The total co-moving length along the line-of-sight of 
20 simulation boxes is \h~^ Gpc, which is about the same size as the observable 21-cm region distributed between redshifts 
10 — 20. The correlation between the boxes can be ignored since the box-size is much larger than the non-linear length scale, 
and the number of neglected modes is small. In Fig. 1, the solid line is the power spectrum of the re-sampled sources. To 
measure the dependence of non-Gaussianity on scale, we compare the results with different scales of experimental noise cut 
off. 

We simply assume the noise to be zero above a cut off and infinity below the cut off scale. This is a reasonable approx- 
imation for a filled aperture experiment, which has good brightness sensitivity, and an exponentially growing noise at small 
scales. Three cut off where chosen at fee = l/iMpc~^, 4/iMpc~^, 16/iMpc~^, which represent the linear, quasi-linear and 
non-linear scales. Three different experimental noise levels are shown as vertical lines in Fig. 1. 

In principle, the convergence map is the variance (or covariance when the filter T'^ has negative value) of the over-density 
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Figure 1. The dimensionloss power spectra of the re-saiiipled dark matter from the 1456"^ N-body simulation in three dimensions are 
given. The solid line is the power spectrum on the 512'^ grids. The re-sampled sources keep the non-linearity and the non-Gaussianity of 
the structures up to fe ~ 30/iMpc~^. Three different experimental noise cut offs are shown with fcc = l/iMpc~^,4/iMpc~^, 16/iMpc~^, 
which represent the linear, quasi-lineaj: and non-linear scales. 



field after a specified filtering process. Sliear is tlie covariance of two maps, since tlie anisotropic filter can not be factored in 
to a perfect square. We need to smootii tiie maps to extract tlie lensing signal with maximum S/N. Tlie window function used 
to smooth the lensed map, which is isotropic in the transverse directions to the hne-of-sight, can be calculated with Eq. (25). 
The gradient of the power spectrum becomes negative at small scales; that comes from the limited resolution of the N-body 
simulation and is unphysical. The experimental noise will put a natural cut off at small scales. 

As mentioned in Section 2.2.2, the reconstruction of k will depend on 2P + AP, in 2-D which is equivalent to the gradient 
of 2-D version of Ain = fc^P2D(fc)/27r. In 3-D, it is more complicated since AP(fc) is not isotropic. The optimal window 
functions have two parts Wi and W2, the choice of which is not unique. One might expect a symmetric decomposition to have 
the best S/N. The optimal filter of k is positive except at a few modes, and can be decomposed in to two equivalent parts 
(one part need to contain a minus sign for those negative value of the filter). In contrast to k, the shear construction needs 
to use the covariance between two different windowed temperature fields, since there is a sin or cos component in the window 
function. The window is a function of the transverse and parallel components of k. 

We can calculate the mean covariance of the two smoothed maps along the redshift axis for each pixel. From Eq. (27) we 
can construct the convergence map. Shear maps are reconstructed in the same way, except different optimal window functions 
are used. The anisotropic part cos20fej^ can be decomposed into cos^fej^ — sinSfej^ and cos^fej^ + sin^fej^. Both windows can 
generate a field with even PDF so that the distribution is less non-Gaussian. This is consistent with the numerical results as 
shown in Fig. 6. Using these two maps, we construct the 71 map with their covariance, as shown in the Eq. (33). Similarly 
we can get the 72 map. 



4 NUMERICAL RESULTS AND DISCUSSION 

Cooray (2004) claims that the variance will not vary considerably and is not a ideal measurement of the lousing signal. Even 
though the k field itself is only a few percent, the integrated effect from the 3-D images will reduce the noise ratio significantly 
to uncover the signal. Zahn & Zaldarriaga (2006) solve the problem from an alternative approach by generalizing the minimum 
variance quadratic estimator (Hu & Okamoto 2002) in CMB lensing to 3-D. 

A related work was done in Benton Metcalf & White (2006), where they also construct quadratic estimators of shear and 
convergence in real space, even though they did not include the correlation between the 2-D slices along the line-of-sight and 
they did not choose the estimator with minimized noise. 




Figure 2. The noise of Icnsing maps from diflferent estimators using experimental noise 1, which cuts off at fee = IhMpc"^. We treat 
the lh~^ Gpc space of gas at z = 10 — 20 as 20 independent sources each is a 50h~^ Mpc box-size cube. Structures at these redshifts are 

similar to those at z = 7 used In' us, though less non-linear. We can expect to see similar non-Gaussianity effects in the reconstruction 
with the Ih^^ Gpc space except that the non-Gaussianity of sources will be smaller. The curves are truncated at V2kc, where the noise 
goes to infinity. The thick solid line is the expected lensing signal. The dotted line is the Icnsing reconstruction noise for a simulated 
Gaussian source with the same power spectrum. The dashed curve is the noise from the N-body simulation using the Gaussian estimator, 
which increases modestly compared to the Gaussian source. It is identical for the optimal k, ■y reconstruction as it is for the deflection 
angle. The thin solid line is noise when shear and convergence are re-weighted by their non-Gaussian variances. 



4.1 Non-Gaussianity 

The dark matter distribution is linear at large scales, and can be treated as Gaussian. In the non-linear scales, when the 
amplitude of density fluctuations is big, the structure becomes highly non-Gaussian. Reference Gaussian sources with identical 
power spectrum to the dark matter are generated. 

We treat the lh~^ Gpc region at 2; = 10 — 20 as 20 independent sources. Structures at these redshifts are similar to 
those at z = 7 used by us, though less non-linear. We can expect to see similar non-Gaussianity effects in the reconstruction 
with the lh~^ Gpc space except that the non-linear scale is smaller. We compare the reconstruction noise with three different 
experimental noise as well as the lensing signal in Fig. 2, 3 and 4. The thick solid line in the middle panel is the lensing power 
spectrum, which is calculated with the Limber integral of the 3-D power spectrum of dark matter using Eq. (4) . We use the 
publicly available code Halofit.f (Smith ct al. 2003) to generate the nonlinear dark matter power spectrum. The code provides 
both their fitting results, and the results using the Peacock-Dodds formula (Peacock & Dodds (1996), PD96 hereafter). The 
'stable clustering' assumption of PD96 breaks down at low redshifts, but is reasonably good at high redshifts where the power 
spectrum is more linear. The halofit code fits the power spectrum at low rcdshift to Virgo and GIF CDM simulations, which 
used the transfer function of Efstathiou et al. (1992). We use a combination of the two: Halofit power spectra are used for 
redshifts lower than z = 3.0, and PD96 power spectra are used for higher redshifts. 

Since the reconstruction noise of k is isotropic, one can always choose the direction of the lensing mode I to be parallel 
with a coordinate axis. In this direction, 71 (Z) = fv(Z),72(i) = 0, and 7e = 71, which simplifies the numerical calculation. The 
optimal combined estimator becomes the sum of k and 7e weighted by their noise. The weights could be the Gaussian k, 7b 
noise, or non-Gaussian noise. We will show that the combined estimator with Gaussian noise weights has the same noise as 
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Figure 3. Same of Fig. 2, but witli cut oflF at tlie quasi-linear scales kc = 4/iMpc~^. The effect of non-Gaussianity of sources is more 
pronounced. We can compare the S/N with a fiducial cosmic shear survey of sources in the same 10 < z < 20 redshift range, which 
reconstructs the Icnsing from the sliapc of galaxies, with a surface density of 14arcmin^-^. To map the Icnsing to the same S/N with 
redshift z ~ 1 sources requires a density of 56 arcmin"^ (Hu & White 2001) with rms ellipticity of 0.4. We see that proposed optical 
lensing surveys are unliliely to outperform 21-cm sources. 

the OQDE for both Gaussian and non-Gaussian sources. Fig. 2, 3, and 4 are results using noise cut offs from experiment 1, 2 
and 3. The curves are truncated at ^/2kc- The non-Gaussianity increased the noise of all estimators. The first cut off falls in 
the linear regime, where the non-Gaussianity only has a modest effect on the noise. The second cut off is at the quasi-linear 
scale. Here the non-Gaussianity increases the noise of the OQDE by about 1 to 2 orders of magnitude. At the highly non-linear 
scales, the non-Gaussian noise is about 3 to 4 magnitude higher than the Gaussian noise, and in fact higher than that for the 
more noisy experiment. 

Our estimators were derived in the limit that k and 7 are constant, and are optimal in that limit. For spatially variable 
lens, we solve for the required normalization factors. In the OQDE, the windows do not depend on the scale of the lens, so 
one might guess the same Ansatz to hold for the (k, 7) estimators. We verify this numerically in Fig. 5. The solid line and 
dotted line is for Gaussian sources and non-Gaussian sources respectively. The differences are less than a few percent, and 
consistent with integration errors from the tabulated power spectrum, and most importantly, independent of scale, as we had 
expected. We do note, that for a finite size survey, the (/t, 7) recover the constant mode, which is lost in the OQDE. Three 
more numbers are recovered. 

The combined estimator with k and 7e weighted by non-Gaussian noise is more optimal than weighted by Gaussian 
noise, therefore has lower noise than the OQDE. In fact, the non-Gaussian noise of 7e is much smaller than k. To investigate 
the origin of this change, we first investigate the cause of the increased noise in non-Gaussian sources for k. This could be 
because either the non-Gaussianity leads to a high kurtosis in k, which boosts the errors; or the non-Gaussianity may lead to 
correlations between modes, resulting in a smaller number of independent modes, and thus a larger error. 

In Fig. 6, the PDF of maps smoothed with the k window are shown. The top, middle and bottom panel show the results 
with experimental noise cut offs 1, 2 and 3. The solid lino is the PDF for maps smoothed with n window (Tf ,T2° in section 
2.2.1). Because the window functions arc almost symmetric, we plot only one PDF. To see the full dynamic range on the 
X-axis, we plot ±|r|^''"' as x-axis, and PDF(|T|^/'*)|T|^®/"' as the y-axis. The integral of the x-axis weighted by the y-axis will 
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Figure 4. Same of Fig. 2, but with cut oflF at the non- linear scales kc = 16/iMpc~^. At the highly non-linear scales, the non-Gaussian 
noise is about 3 to 4 magnitude higher than the Gaussian noise. The combined re- weighted estimator (NG-NG k + 7e ) has noise about 
half an order of magnitude lower than the OQDE. 



give (r''), which is basically a estimation of the point-wise non-Gaussian reconstruction noise. Here VT>¥{\T\^'^) is the PDF 
of \T\^'^. To compare with a Gaussian distribution, dotted lines are also plotted. The contributions to the (T*) in experiment 
1 mainly come from small fluctuation regions. In experiment 2, the large outliers play a more important role but one can still 
expect the curve to converge. In experiment 3, most contributions come from rare regions with high fluctuations. Gaution 
should be exercised in the interpretation of the most non-linear scales, since a larger number of source samples may result in 
a different error. It is clear, however, that the noise has increased dramatically. 

The kurtosis of k is {{T^f)/{{T^ff - 3, and an analogous quantity can be defined by {{T^T^f)/{{{T^f){{T^f)) - 1 
for shear. Tf ~ T2 , and is uncorrelated with 7^7 . The noise of k and 7 is determined by both kurtosis and number of 
independent cells. For experimental noise 1, the kurtosis of T" and are 1.2 and 0.29 respectively. The effectively independent 
cube cells for k and 7 have side length 4.8/i~^ Mpc and 4.6/i~^ Mpc respectively. The corresponding Gaussian sources with 
the same cut off have effective cell size 3.0/i~^ Mpc and 3.5/i~^ Mpc. For experimental noise 2, the kurtosis of and are 
18 and 5.7 respectively. The effective cell size for n and 7 are 1.8/i^^ Mpc and 1.5/i~^ Mpc respectively. The corresponding 
Gaussian sources with the same cut off have effective cell size l.Qh~^ Mpc and l.lh~^ Mpc. For experimental noise 3, the 
kurtosis for T'^ and T'*' are 1.6 x 10^ and 3.5 x lO'^ respectively. The effective cell size for k and 7 are 540/i~^ Kpc and 
310/i^^ Kpc respectively. The corresponding Gaussian sources with the same cut off have effective cell size 240/1" ^ Kpc and 
290/i~^ Kpc. We conclude that the shear measurements have lower non-Gaussian noise both because of a smaller point-wise 
kurtosis, and less correlation between modes. 

We will see later that experiment 2 has the largest S/N, which is larger than unity for I < 6000. We can compare the 
S/N with cosmic shear surveys, which reconstruct lensing from the shape of galaxies. The noise can be estimated by {■y^)/ne« 
(Hoekstra et al. 2006; Hu & White 2001), where we use rms intrinsic ellipticity, and UcS is the effective 

number density of galaxies. We plot the shear noise from a survey of sources in the same redshift range 10 < z < 20 in Fig. 3, 
with a surface density of 14arcmin~2. For more realistic source redshifts z ~ 1 in proposed optical surveys (Hu & White 
2001), this corresponds to a surface density of SGarcmin"^. in the CFHTLS wide survey the source galaxies are distributed 
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Figure 5. The comparison of reconstruction noise from the combined (ft, 7) estimator and the OQDE. While the optimality is only 
proved at low I, we find them equally optimal for Gaussian sources at all scales. The scatter is consistent with numerical integration 
errors from the tabulated power spectrum. 



at redshifts lower than 3, and their effective number density is ~ 12 galaxies arcmin"^ (Hoekstra et al. 2006). This noise 
is larger still. Even though non-Gaussian 21-cm lensing saturates lensing reconstruction, it still measures more modes than 
current proposed optical surveys. 

In Fig. 7, we show the reconstruction noise at two different / versus various oxpcrimcutal noise cut off fee. The top panel is 
for the fundamental mode in the box, Zi = 2n/L = 783, and the bottom panel is for I2 = 6I1 = 4715. As shown in the plot, it 
is clear that the noise of a Gaussian source decreases as fee increases, because of the increasing number of independent modes. 
The dotted lines are a least squares fitting power law Nok~^ to the Gaussian noises, and A'o — 3.1 x lO^'^, 1.3 x 10~^ for top 
and bottom panels respectively. This comes from counting the number of available source modes. The dashed lines connect the 
non-Gaussian noises of the OQDE. The triangles are the reconstruction noise for the combination estimator, which is equal 
to the OQDE at larger scale kc and about half an order of magnitude lower at smaller scales of kc- From this plot, wc can 
see that experiment with lower noise does not necessarily decrease the reconstruction noise of the OQDE for non-Gaussian 
sources. And the experimental noise has a limit around the quasi-linear scale where the OQDE achieves its best S/N. The 
S/N achieves its maximum around kf^' « 4/iMpc~^. This cut off with maximum S/N varies only slowly with /. 

If one wants to estimate the effective number of available lensing modes, we can derive an effective cut off of a Gaussian 
field which gives the same S/N as the optimal non-Gaussian sources estimator. This is « 2/iMpc~^, where the power 
spectrum of source is ^ 0.2. The size of the effectively independent cells is 2.07i ^ Mpc. A simple Gaussian noise estimate 
counts all modes up to A^(fc) < 0.2, which is perhaps surprisingly low. 

For our noise estimates, we stacked simulations all at redshift z = 7. While the angular diameter distance does not change 
much to ~ 20, the structure does evolve. We do not have access to the higher redshift outputs to test this effect, but one 
would expect a smaller non-linear scale to result in a smaller reconstruction noise. 
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Figure 6. The PDF of maps smoothed with k window are shown. The top, middle and bottom panel show the results with experimental 
noise cut offs 1, 2 and 3. The solid line is the PDF for maps smoothed with the k window (T"^ in section 2.2.1). To see the full dynamic 

range on the x-axis, wc plot the curve with ztlT]-'/^ as x-axis, and PDF(|T|-'-/'*)|T|-'-'''''^ as the y-axis. The integral of the x-axis weighted by 
the y-axis will give (T^), which is basically a estimation of the reconstruction noise. The error bars are estimated from the 20 simulations. 
To compare with a Gaussian distribution, dotted lines are also plotted. The contributions to the (T^) in experiment 1 mainly come from 
small fluctuation regions. In experiment 2, the large outliers play a more important role but one can still expect the curve to converge. 
In experiment 3, most contributions come from rare regions with high fluctuations. Caution should be exercised in the interpretation of 
the most non-linear scales, since a larger number of source samples may result in a different error. It is clear, however, that the noise has 
increased dramatically. 



4.2 Future directions 

A possible way to find the optimal window functions for non-Gaussian sources is to divide the window into A'^ frequency 
bins Wi(fci, fc2, fejv), and apply a numerical variation to those bins. The noise can be measured numerically by applying 
the estimator to the simulated sources. The process of searching for a optimal filter is equivalent to look for a minimum of 
reconstruction noise in A^ dimensional space fei, fe2, fciv. One could use a Newton-Raphson method to do this. In this paper 
we only considered the class of windows which are identical to the optimal Gaussian estimators with a hard cut off, as well 
as two weightings for shear and convergence. 

One can also try to Gaussianizc the sources by modifying the PDF of all the sources to be Gaussian. The physical 
explanation and details of Gaussianization can be found in Weinberg (1992). The basic idea is that every pixel should preserve 
its rank in the whole field during the Gaussianization process. During structure formation, the non-linear evolution at small 
scales shotild not destroy most of the information on the peaks and dips of the linear field. However, this Gaussianization 
process will change the power spectrum of sources, and the reconstructed lensing field will be biased. This is not a linear 
process, and the variation of power spectrum does not have analytical solution, and can only be measured numerically with 
simulated sources. 

Recently it has been proposed that one could economically achieve brightness mapping of 21-cm emission at lower redshifts 
(Chang et al. 2007), potentially even with existing telescopes. If individual galaxies are not resolved, one can again ask the 
question of how one could reconstruct a lensing signal. This is very similar to the problem studied in this paper. 
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Figure 7. The reconstruction noise versus the cut off in the experimental noise. The top panel is for li = 27r/L = 783, and the bottom 
panel is for I2 = 6h = 4715. The noise of Gaussian sources decreases as fee increases, because of the increasing number of independent 

modes. Tlic dotted lines arc a least squares fitting power law N^k^^ to the Gaussian noises, and A^o = 3.1 X 10~^, 1.3 X 10^-'^ for top and 
bottom panels respectively. The dashed lines connect the non-Gaussian noise of the OQDE. The triangles are the reconstruction noise 
for the combined estimator, which is equal to the OQDE at larger scale fcc and about half an order of magnitude lower at large kc- The 
noise of the non-Gaussian sources changes slowly and saturates or even increases at small scales. 



5 CONCLUSION 

In this paper, we developed the maximum likelihood estimator for the large-scale structure from the 21-cm emission of the 
neutral gas before the epoch of re-ionization. The convergence and shears can be constructed independently. To test the effects 
of non-Gaussianity, we applied our estimators to simulated data. The sources wore generated by N-body simulations, because 
gas is expected to trace the total mass distribution. To investigate the influence of non-Gaussianity, we also use Gaussian 
sources which have the same power spectrum as the simulated sources. We applied our estimator and the OQDE on both the 
Gaussian and non-Gaussian sources. Though our estimators are derived in the simplified case of a constant convergence, the 
noise of our combined estimator of convergence and shear are the same as the OQDE for Gaussian sources. For a finite survey 
area, three extra constant modes can be recovered. 

The non-Gaussian nature of the source can increase the error bar by orders of magnitude, depending on the experimental 
cut off scale. Shear construction is affected less by non-Gaussianity than the convergence field, and the combined estimator 
with non-Gaussian noise weights is a better choice than reconstructing with the OQDE. S/N can not be boosted infinitely 
by reducing the experimental noise, and achieves its maximum for a cut off around k^'^ « 4/iMpc~^. Below that scale 
the S/N start to saturate or even decrease. The maximum S/N for non-Gaussian sources is equal to Gaussian sources with 
« 2/iMpc~^, where the power spectrum of source is « 0.2 and the side length of the effectively independent cells is 
2.0/1"^ Mpc. The maximum S/N is greater than unity for / < 6000, which makes 21-cm lensing very competitive compared to 
optical approaches. 
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APPENDIX A: NORMALIZATION AND NOISE OF THE ESTIMATOR 

In the end of section 4, the numerical results of the noise of the estimators are shown. Here we will develop the analytical 
expression for 

S«(fe±) = 6«(fex)[/t(fe±)+n(fex)]. (Al) 

For shear, a similar relationship holds even though b and n are not isotropic. 

When K is spatially variable, 

ft{x) = Tb(x^- D{xi_),x\\) = Tt{x±,x^\) -V±Tt{x±,x^\) ■ D{x±) , (A2) 

where D{x±) = d{x±)xizs), and d{x±) is the deflection angle. Therefore k = Vx • D- 
Fourier transforming Eq. (27), 

= J d\±E4xi_)e-"'^-'^^ = ^ J d'xf:^{x)f:^ix)e'"'^-^^ - {2nf6^''{k^)V^ . (A3) 

T = Tb + n, and noise is uncorrelated with the signal. The product in real space can be represented as a convolution in Fourier 
space 

I d?xe-"'—n^{x)n^{x) = I ^T:^(fe;,fef|)f:,(fcx -fe'x,-fcf|). (A4) 

fi,(fe)= J d\e-""''n{x±-D{x±),xii)=Tt,{k)- J d\e-"'""V±n{x±,xii)-D{x±), (A5) 
and the lensing introduced term can be further simplified as 

J d'xe-"' ''V±n{x±,x^^) ■ £»(xx) = y d'a;e-"= "Tt,(xx,a;||)(ifex - Vx) • £»(xx) 
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The quadratic term in Eq. (A3) can be written as 



n{k± - fe'x,fe||)(ifex - Vx) ■ i3(fe'x) . (A6) 



J d\e-"'^-^^f:^{x^,x\Of:,{x±,x\\) = j A^l^f(feY,^)W'2''(fcx-fcl,-fci|)[7^b(fc'x,fe'|)T,(fcx-fc'x,-M|) 

-T6(fex-fel,-fef|) j i^T6(fe'x-/cl,fef|)(ife'x-Vx)-£>(fel) 

-T6(fcl,fc||) J ^r,(fex - fe'x - kr,-k\) (i(fex - fex) - Vx) • -D(feT)] 
+Noise . (A7) 

Using the relationship that 

(Tb (fe'x, fcf|)T6 (fex - fe'x, -fey)) = (27r)3,53°(fex,0)P3D(fex,fe||) , (A8) 

we found that the expectation value of the first terms and the noise term in Eq. (A7) can cancel the last term in Eq. (A3). 

Note (5(0) = limAfc^o(Afc)-i ~ {L/2n), and W^{k± - fe'x, -fe||) ~ W^2'(fe'x, fey) since <5^°(fex) is nonzero only when fex = 0. 
Similarly, the last two terms can be simplified. The normalization factor 

b4k±) = ^ J ^^l^r(fe'x,fe||)W'2''(fex-feY,-feI|)[(fex-fel)-fexP3D(fex-fel, -fell) + fe'x- fex P3D(fe'x, fey)]. (A9) 

Similarly, replacing W^, W5 by W^\W^^ {WJ\W^^), and fei by fei cos26lfe^ (fei sin26lfe^), we find the normalization factor 
for 71 (72). 

The noise of the estimator can be calculated in the absence of lensing: (|«;(fex)P) = («:(fex)«;*(fcx))- Since {jK(fex)P) = 
(27r)^52°(0)iV«(fex) and <5^°(0) = limAfc^o(Afe)-2 ~ (-C'/27r)^ Wick's theorem gives 

^"^^"-^ = KfeIpL y (^{■f'»(fc^-fcx,-fey)P3D(fc'x,fci|)rf(fcx-feY,-fe||)H/^(fe'x,fey)]' 

+ P3D(fcx -feY,-M|)f'3D(fc'x,fc'|)J^"(fc± -fc'x,-fc'|)J^"(fcY,fe'|)}- (AlO) 

The first term is the convolution of P3D(fe)W^i''(fe)^ and P3D(fe)W2''(fc)^, and the second term is the convolution of P3D(fc).^''(fc) 
with itself. The dimensionless quantity fex-^«(fex)/(27r) is equivalent to PCi/(27r) in other literature. 
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